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ABSTRACT 



This work is directed toward improving an existing general 
circulation model (Haney 1974) through the inclusion of time 
dependent surface boundary conditions and a wind mixing 
parameterization . 

Climatological atmospheric values, functions of time, 
provided the model the needed parameters in order to calculate 
solar insolation, longwave back radiation from the ocean, 
sensible heat exchange and latent heat exchange, The mixing 
parameterization involved the use of the Monin-Obukov length 
scale as a measure of the mixed-layer depth. The first law 
of thermodynamics was used to calculate the form of eddy flux 
above the mixed layer depth, 

A one-dimensional model using two different types of 
boundary conditions was used to test various parameters and 
ideas. The boundary condition similar to that used in the 
baroclinic model produced results through an interacting 
system that appeared similar to observation at Ocean Station 
"N," 

The three-dimensional model was integrated through 10 
years of time. Results showed a lack of a paramanent 
thermocline due to a strong vertical diffusion. A lack of 
behavior attributable to mixing was noticed, indicating the 
climatological winds were too weak for mixing. 
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I . INTRODUCTION 



The possibilities and ramifications of accurate, long- 
range weather forecasts are extremely important to the United 
States and the world, There have been annual changes in 
global atmospheric weather patterns documented for some 
time reflected by changes in temperature, rainfall and wind 
patterns. Similarly, the ocean experiences variations in 
circulation and temperature patterns that influence fluxes 
across the air-sea interface and thus affect the atmospheric 
weather patterns, 

Temperature variations appear on vertical and horizontal 
scales of hundreds of meters in depth, millions of square 
kilometers in area, and contain temperature variances of 1°C 
to 2°C above or below normal. As the fluxes across the air- 
sea interface of such an area represent enormous amounts of 
thermal energy, the presence or absence of such temperature 
anomalies may well influence the heat and water exchange on 
a scale sufficient to modify existing weather patterns. It 
is not certain whether the thermally anomalous patches can 
initiate shifts in climate or if they are a passive response 
to the atmosphere. 

In the North Pacific Experiment (NORPAX), scientists are 
attempting to learn about these patches. The Main goals of 
NORPAX are to learn how the patches are formed and what their 
role might be in the exchange of heat and vapor to the 
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atmosphere. In the long run, scientists hope to be able to 
monitor ocean temperature variations by remote sensing and 
then to use that data to predict some features of weather 
patterns and climatic shifts months and years in advance. 

Although reliable long range prediction is still many 
years away, J. Namias of Scripps Institute of Oceanography 
and other NORPAX scientists have discovered some interesting 
correlations. In many years when very intense sea surface 
temperature anomalies were present in the North Pacific, 
there were also distinct weather patterns present over the 
United States (Namias, 1969). Although the results are not 
completely conclusive, they are suggestive of a link between 
the Pacific temperature anomalies and the climate downwind 
over the United States. 

The NORPAX work is concentrated in several areas, however, 
only the theoretical/numerical area is pertinent to this 
thesis. Numerical models attempt to simulate the ocean 
dynamics and numerically integrate in time to form predictions. 
Since NORPAX deals with temperature anomalies on the order of 
1°C - 2°C, accuracy in the description of upper ocean dynamics 
is essential, Thus all advection, diffusion, convection, and 
wind mixing should be as accurately portrayed as possible. 

The numerical model on which this work was based is that 
of Haney (1974). This was a 6-level model of a baroclinic 
ocean with flat bottom and rectangular coastline from 51,25°S 
latitude to 48.75°N latitude. The primary goal of this work 
was to make improvements in his model so it could be used to 
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predict the formation and evolution of the large-scale, sea- 
surface temperature anomalies described above. The improve- 
ments made were of two general types: improvements in 

resolution and improvements in parameterizing the physical 
processes . 

In order to improve accuracy in the model, especially in 
the upper layers, Haney's model was expanded to 10 layers. 

The layers were made thinner in the upper ocean for resolution 
considerations, In addition, the geographical domain was 
reduced to include that part of the North Pacific Ocean 
bounded by the equator and 65°N, Of course decreasing grid 
sizes alone could not suffice to enhance accuracy. The 
dynamics of the model had to be modified to provide appropriate 
parameterizations of the pertinent physical processes, 
Improvements were therefore made in calculating the surface 
fluxes and in parameterizing the effects of wind mixing. 

Haney's 1974 model used a simple linearized equation to 
calculate the heat fluxes across the air-sea interface; 

FLUX = K (T a - T gfc ) (1) 



where T^ was an equilibrium air temperature and was 

the predicted sea-surface temperature. K was a coefficient 
estimating the combined effects of sensible heat flux, upward 
long-wave radiation, latent heat release and downward short- 
wave solar radiation. In the improved version of Haney's 
model, each of the above surface fluxes was individually 
computed using prescribed atmospheric data, The atmospheric 
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data included incoming solar radiation, cloud coverage, 
surface air temperature, vapor pressure, and both north- 
south and east-west components of the surface geostrophic 
wind. In the actual time integration of the model the fluxes 
were at first computed using annual mean atmospheric param- 
eters. After initial shocks in the model decreased, the 
first two annual harmonics of the atmospheric parameters 
were introduced to give time-dependent boundary conditions. 

Haney's model provided vertical sub-grid scale heat 
fluxes through diffusion and convective adjustment. Vertical 
eddy diffusion of heat was accomplished with a constant 
diffusion coefficient whose magnitude was appropriate for 
the deep ocean. This process was retained along with the 
convective adjustment, but a new wind-mixing parameterization 
was introduced. According to recent mixed layer theories, 
(Denman, 1973), above a certain mixed-layer depth, wind 
generated turbulence and surface convections should both 
serve to mix heat downward across the interface. With finer 
vertical resolution in the improved model, vertical eddy 
diffusion using a constant coefficient was not adequate. 

The surface stress used in Haney's (1974) model was an 
idealized one based on the mean zonal wind stress for all 
oceans combined (Hellerman, 1968). In the improved model, 
data from the North Pacific Ocean was used. As with the 
improved heat flux scheme, annual mean stresses were used at 
first until the initial shocks decreased. Then the first 
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two annual harmonics of the winds were added to give time 
dependency. Thus, the stress supporting both the circulation 
and the mixing were made time dependent. 
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II. 



THE NUMERICAL MODEL 



A. RESOLUTION AND COMPUTATIONAL EFFICIENCY 

As stated above, the basic model used in this work was 
that of Haney (1974). Many improvements were made in the 
model to effect the desired modifications, "The model is 
based on the hydrostatic and Boussinesq approximations. This 
means that the ocean is assumed incompressible, and the 
density is replaced by a constant everywhere except in the 
hydrostatic equation where it is multiplied by the acceler- 
ation of gravity, A major simplification is the neglect of 
salinity and the assumption that the density is a linear 
function of temperature alone,"'*' 

The reader is referred to Haney's 1974 publication for 
details but in the interest of completeness the following 
basic equations are reproduced. The equations of motion; 



du 

dt 



-1 



3P 



p Q a cos <f>3A 



+ fv + 



u 



a cos <J) 
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dv = In _ f 

dt p a 3<f> 
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u 

a cos <j> 
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sin <J> ■+ A m ^ 2 v 



+ k 



3 2 v 
3z 2 ' 



(3) 



Haney, Robert L. , "A Numerical Study of the Response 
of an Idealized Ocean to Large-Scale Heat and Momentum Flux," 
Journal of Physical Oceanography , v. 4, No. 2, p. 147, 1974. 
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The hydrostatic equation; 



3p 

Iz = -PS 



( 4 ) 



The mass continuity equation; 



3w 

3z a cos 



i r 3u 3 . , ri ,, 

sf* h (v cos = 0 



(5) 



The first law of thermodynamics; 



f - A H V2T + k + S c< T > 



( 6 ) 



The equation of state; 



p = p o [1 - a(T - T o )] , 



(7) 



The symbols are listed in Table (1) for clarity, 



TABLE 1 



SYMBOL 



DEFINITION 



t 

X 

<*> 

z 

ft 

f 

a 



g 

u 

V 

w 

T 



P 

P 

a 



time 

longitude 

latitude 

depth 

angular speed of the earth's rotation 
Coriolis parameter ( = 2 sinft) 
radius of the earth 
acceleration of gravity 
eastward velocity component 
northward velocity component 
vertical velocity component 
temperature 

constant reference temperature 
density 

density of water at reference temperature 
coefficient of thermal expansion, 
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Several improvements were introduce to improve the 
resolution of the model. Since the work was expressly 
designed to improve accuracy in the sea-surface temperature, 
the vertical resolution had to be adjusted. The model was 
vertically expanded to ten levels fixed at 10, 30, 60, 100, 
150, 225, 350, 700, 1500, and 3000 meters. These locations 
provided greater precision in the top-most layers of the 
ocean. This was in consonance with the effort to better 
parameterize the effects of wind and vertical motions which 
are so important in changing the temperature in the upper 
layers of the ocean, The total depth of the ocean remained 
constant at 4000 meters. 

The horizontal gridding was also altered. In the original 
model, the latitude-longitude grid intervals were 2.5° and 
3.0^ respectively. In the present model, the latitude- 
longitude grid intervals were 2.0 < ’ > and 2.8° respectively. 

Since the actual location of the domain was altered, the 
lateral boundary conditions were altered at the southern 
boundai'y (equator) where a free slip (no friction) condition 
was imposed. 

Many schemes for improving the computational efficiency 
of the model were attempted; those having favorable trade- 
offs between computer space and computer time were retained. 

By far the most successful ploy involved a very fast Poisson 
equation solver designed by Roland Sweet (1972). The model 
used a streamfunction ¥ to predict the vertically averaged 
mean currents (u, v), The streamfunction tendency was 
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previously obtained by solving the Poisson-type equation 
every time step by over-relaxation techniques, Even though 
an optimum over-relaxation factor had been used, the use of 
the direct solver produced a savings in computer time of 
approximately 30%, This figure, however, is variable depend- 
ing on the portion of the program dealing with the solution 
of the Poisson equation. If a model spends a large percentage 
of its time in solving the Poisson equation, a much greater 
percent of savings can be realized, A barotropic model or a 
baroclinic model having fewer vertical levels would be such 
an example. 

An attempt at eliminating repeated multiplications every 
time step was made. Because of the spherical grid, all 
transports of heat and momentum had to be weighted by their 
appropriate area or volume. The various combinations of grid 
sizes (e.g. Ax X Ay) were stored and called from core when 
used vice their repeated multiplication. It was found that 
this ploy made insignificant savings in time while expanding 
computer space greatly. It was thought the time spent dealing 
with the indices associated with the storage was comparable 
to the multiplications, This may not be the case on a 
different computer. All work discussed herein was performed 
on the IBM 360. 

Another attempt at economizing time was to eliminate the 
need to switch the fields associated with one time (say t^) 
and the previous time (t^ -At) before calculating the 
friction terms. In the original model, this switching was 
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done every time step except every tenth step which employed 
a Matsuno scheme, This was altered so that the fields were 
switched only at the infrequent Matsuno steps. Again the 
time savings was small in that only a small portion of time 
was spent doing this operation, However the scheme did not 
create a need for storage expansion and did seem more logical 
in its notation and was therefore retained, 

B. TIME DEPENDENT BOUNDARY CONDITIONS 

In addition to improvements in resolution and efficiency, 
improvements were made in representing the physical processes. 
These changes were in two main areas, the first of which was 
in the calculation of the fluxes of heat and momentum at the 
sea surface. 

As shown in Eq ,( 1 ), Haney ' s 1974 model employed a simple 
linearized version of heat flux. In an earlier paper, Haney 
(1971), formulated the flux by separately calculating the 
incoming solar radiation (Qj), the net upward longwave 
radiation (Qg), the latent heat exchange (Qg)> and sensible 
heat exchange (Qg). The net downward heat flux Q was then; 

Q = Q I - ( q b + q e + Q h ) • ( 8 ) 

r* 

The present model employed this method of calculating the 
heat flux. Back radiation was calculated from 

Q b (T 0 ) = Q* o T o - (9) 

where a was the Stef an-Boltzman constant, T q was the 
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sea^surface temperature calculated by the model, and 

Q* = .985 (.39 - ,05 e *) ( 1 - ,6C 2 ) (10) 

a 

where e a was the atmospheric vapor pressure at 10 meters 

and C was the fraction of cloud coverage. This dependence 

* 

of Q on vapor pressure and cloud coverage parameterized 
the reduction of upward longwave radiation by water vapor 
and clouds in the atmosphere. 

The latent heat flux (CL,) from the ocean was 




where p„ was the density of the air, CL. was the drag coeffi- 

ct ]J 

cient, | V | the wind speed at anenometer level, L the latent 
heat of vaporization and the atmospheric pressure. The 
vapor pressures e s (T 0 ) and e a were the saturation vapor 
pressure at the sea surface temperature and the atmospheric 
vapor pressure at 10 meters, respectively. 

Sensible heat exchange (Q„) across the interface was 

xi 

given by 



«H * 



P C. 
a D 



|V| c 



(T - T ) 
v o a y 



( 12 ) 



where was the specific heat of air at constant pressure 
and T^ was atmospheric temperature at 10 meters. 

Incoming solar radiation (Qj) was calculated with 



Qj = S Q (.74 - ,6C) . 



( 13 ) 



This expression parameterized the albedo of a cloudless 



16 



atmosphere due to air molecules, dust and water vapor scatter- 
ing plus the albedo of clouds, S q was the solar insolation 
incident on a horizontal surface at the top of the atmosphere. 
Parameters not discussed in Equations (8) - (13) are given 
in Table (1), The above heat flux components were calculated 
at each time step from the sea surface temperature predicted 
in the model and the atmospheric parameters which were 
specified on a time-dependent basis, 

In addition to the thermal forcing, the ocean circulation 
was also forced by surface winds given on a time-dependent 
basis. Assuming an average geostrophic inflow angle of 10° 
and a frictional reduction of velocity by 10%, one could 
derive the surface stress from the surface geostrophic winds, 
if u and v are the eastward and northward surface winds, 
respectively, after inflow and frictional influences are 
considered, then 

| V | = (u 2 + v 2 )^ 

and the stresses 

T x = p a C D l y l u (14a) 

T y = p a C D v (14b) 

The atmospheric data representing monthly climatology was 
obtained through various sources. Solar insolation values 
were obtained from the Smithsonian Meteorological Tables 
(List, 1963, Table 132), Monthly values of the cloud coverage 
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were obtained from the satellite based cloud atlas of Miller 
(1971), Monthly values of the surface (10 meters) air 
temperature, vapor pressure and zonal and meridional geo- 
strophic wind components were all obtained from an NCAR data 
tape of the Northern Hemisphere Climatological Atlases 
of Jenne ejt. al, (1974). 

In order for the ocean model to simulate a realistic 
monthly normal climatology, the atmospheric wind and thermal 
forcing which drives the circulation must vary in a continuous 
manner over the annual period, To provide this continuity 
in the atmospheric forcing, the annual mean and the amplitude 
and phase of the first two harmonics were calculated from the 
above climatological data. The data was averaged along each 
grid latitudinal line with the exception of v, the northward 
component of the surface geostrophic wind, which was averaged 
along each meridion. 

Thus, the atmospheric forcing was completely specified 
on a continuous time dependent basis. In operation, the 
atmospheric variables were employed in two stages. The first 
stage used only the annual mean data. Starting from conditions 
of rest and temperature a function of <J> and z only, the 
integration was continued for about 5 years until the initial 
oscillation "settled" down, Initial shock was severe due 
mainly to large shear currents and boundary effects. The 
model was then subjected to the time dependency in the 
atmospheric boundary conditions. It was expected that 
temperature change, in the upper few hundred meters, due to 
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diffusion from the surface would occur within 10 years, 

Also, the initialization was thought to be a good estimate of 
actual conditions. Thus, since the depth of interest was 
the upper ocean, and the seasonal variation would be dominant 
there, a further integration was not carried out, 

C, THE MIXING PARAMETERIZATION 

The most important change in the model was a parameter- 
ization of the effects of wind generated mixing based upon 
the mixed layer theories of Kraus and Turner (1967) and 
Denman (1973), The change was two-fold in that both the 
depth to which mixing was to take place and the form of the 
mixing needed to be parameterized. 

The parameterization of vertical eddy fluxes in the mixed 
layer involved the critical assumption that the Monin-Obukov 
(M-0) length scale determined the depth in the ocean to 
which mixing occurred, This depth is generally understood 
to be the depth at which mixing due to buoyancy forces 
(convection) is of the same magnitude as the mixing which is 
mechanically driven. At a depth less than the M-0 length, 
mechanical mixing would dominate whereas at a depth greater 
than the M-0 length, mixing would be predominately due to 
buoyancy forcing, 

From Phillips (1966), the equation for the M-0 length, 
with buoyancy fluxes due to salinity neglected, was 

h = w*/(k g a Q) (15) 

where 
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( 15 ) 




and k was the von Karman constant, 

The net downward heat flux Q could have been extremely 
small in some cases. This would have led to an inordinately 
large value of h, In the model however, this mixed layer 
depth was not allowed to extend below the top of the model's 
bottom layer, A more realistic procedure would have been 
to constrain the depth to less than w + /f where f was the 
Coriolis parameter. Both constraints were attempted with no 
ill effects (see results). In operation 

h = min (h , ) 

v mo' f } 

where h was the M-0 length, This gave the depth to 
mo 

which surface originated mixing extended. 

The question then was what form the mixing would take. 

To determine this, the first law of thermodynamics was 
written : 

|| + V-(VT) + f^(wT) = A V 2 T + || - || + 6 c (T) (17) 

where 

F = - k + (w 7 ! 71 ") 

d Z 

This was similar to Haney's (1974) equation with the 
addition of the (w'T') term representing the vertical mixing 
of heat due to eddy fluxes. Although in the model, T and V 
were not necessarily independent of z above h, this was 
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assumed for the purpose of calculating the wind mixing term, 

Also for this purpose the convective adjustment term 6(T) 

c 

was assumed to be negligible. Thus, equation (17) could be 
written 



9T 

at 



+ 



V* (VT) 



9z 



(wT) - AV 2 T - k 



a 2 T 

az 



- 6 (T) = Tp- 
c ' az 



a (w'T’ ) 

az 



and by the above assumptions the left hand side was indepen- 
dent of z. This result is comparable to those of the one 
dimensional mixed layer theories in which advection, diffusion 
and convective adjustment are neglected altogether, and T is 
assumed independent of z in the mixed layer. Thus 



as 3(w'T' ) _ , 

az az D 



(b = constant) 



and integrating to an arbitrary z 



S(z) - (w'T') = a + bz (linear in z). 

z 



Applying the boundary conditions a z=0 and z=-h, the constants 
a and b were determined: 



for z=0 
z=-h 



a = S Q -(w'T* ) Q 
a - bh = S_ h -(¥ r T r )_ h 



or 

b - E a - <S -h '< i » T T T >-h> 



Substituting into the integrated thermodynamic equation 
one obtains 
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(w'T' ) = (w'T' ) + z/, |'(w 7 'T T ) 
v 7 z v ' o ' h L v J o 



(w’T ' )_ h ]+ S(z) 




( 18 ) 



The expression for (w'T 1 ) came from the vertical integral 
of Eq , ( 18 ) and use of the steady state turbulent kinetic 
energy equation 




(w'T’ ) 

z 



dz 



D 



G 



(19) 



where D and G are dissipation and generation of turbulent 
kinetic energy respectively, Integrating Eq,(18) from (-h) 
to 0, and using Eq.(19), the resultant equation was solved 
for (w'T’ )_ h : 



1 2 

(w'T')_h = -2/ h (G - D) - (w'T’) o + 

( 20 ) 

The terms on the right hand side were the contributions to 
thermal flux through the mixed layer at -h. Term 1 represented 
the contribution of wind mixing, term 2 that of convective 
mixing and term 3 the stabilizing effect due to absorption 
of solar radiation. The effects of solar heating were 
interpreted as follows by Han (private communication), If 
the ocean were completely opaque to solar radiation the result 
would enhance stability. If the flux of solar radiation had 
a linear profile with depth, there would have been no effect 
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on stability. This was due to a uniform heating throughout 
the layer (see Fig, *1 ) , If the flux took on an exponential 
profile with depth, the effect would have been one enhancing 
thermal stability. 




Figure 1. The stabilizing effect of solar radiation. 



Paulson (1974) recently found that the profile of solar 
flux in the ocean could be satisfactorily represented 
exponentially. On a scale consistent with this model, he 
found 

S(z) = S Q e (6 ’ z) 
where z < 0 and 3' = (l/19)m 1 . 
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The generation of turbulent kinetic energy similar to 



Denman (1973) took the form 



G = 



3w 



3 

*/ag 



where 3 ~ 1.0, The dissipation D was assumed zero. 

Thus, the defining equations for (w' T' ) , and (w' T' ) 

n z 

were formulated: 



(w'T')_ h = ($w| ) - ( w ' T 1 ) + S + s_ h 

1 ag 

i^ S o - S -h> 



( 21 ) 



and, thus the general expression for the heat flux may be 
written : 



(w'T') z = (w’T’) o + z/ h [(w'T') o - (w’T’)_ h ] 

+ S(z) - S Q - z/ h (S Q _ S _ h ) (22) 

If the value of ( w ' T 1 ) ^ calculated from Eq , (21) was 
positive, indicating upward heat flux at -h, then (w'T* )_^ 
was set to zero. This criteria was set upon (w'T')_ h 
because the interface below the mixed layer is essentially 
thermally stratified. Wind mixing would tend to reduce the 
stratification, never giving an upward flux of heat, 



24 



Ill , RESULTS 



A, TYPE I BOUNDARY CONDITIONS 

In order to test various aspects of the mixing parameter- 
ization before applying it to the three-dimensional, 
baroclinic model, a one-dimensional model was devised, In 
this model, horizontal variations were neglected and there- 
fore no vertical motions due to continuity equation consider- 
ations were present, The consequence was a vertical "pole" 
of ocean in which temperature was affected by flux across 
the air-sea interface, convection, wind mixing and vertical 
diffusion . 

The boundary conditions on the interface flux were 
prescribed in two ways, In TYPE I, a time-dependent heat 
flux was specified (Fig, 2) guaranteeing that the integral 
of interface flux over the annual cycle was zero. This 
ensured that the model would neither heat nor cool excessively. 
In TYPE II, a time-dependent air temperature was prescribed 
from which the heat flux was determined, The condition did 
not ensure a zero annual average but ensured that the "pole" 
ocean would seek an equilibrium heat content and distribution 
that would eventually repeat itself over the annual cycle. 

Neglecting horizontal variations and motions, the first 
law of thermodynamics was written 
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where 



F = -k + (w ' T ' ) , 

d Z 

Thus the change in temperature was due to four physical 
processes: diffusion, convection, eddy flux due to wind 

mixing, and surface eddy flux. 

Many tests were made using both types of boundary 
conditions. In several instances one physical process was 
omitted to test its significance. Other tests included 
modifications of constants employed and variables prescribed, 
After evaluation, a standard run was selected on the basis 
of comparisons with other runs and observations. Figures 3a 
and 3b represent the standard TYPE I run, and show temperature 
as a function of time and depth. The ocean was initially 
isothermal at 20°C down to a depth of 4000 meters. Since 
the changes in temperature were surface oriented, and since 
there were no vertical motions, there existed no permanent 
thermocline. Thus the depth of interest was somewhat shallow 
and the graph extends to only 280 meters. As in all runs 
made, the graph represents the twentieth year of integration, 
time increasing left to right. 

As seen in Fig, 3a there exists a basic asymmetry 
between the spring and fall periods. In spring, surface heat 
flux is downward tending to warm the ocean. However, only 
diffusion and wind mixing transport the heat down. Since 
h is extremely small, the heating is confined to the upper 
layers. In fall heat flux is upward tending to cool the 
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surface layers. However, in this case convection causes the 
entire column to be cooled forming the abrupt isothermal 
layers. The temperature variation lags the boundary condition 
heat flux input by approximately 2,0 months, The maximum 
heat flux occurs July 1 and the maximum sea-surface temperature 
occurs about August 25, Figure 3b is a series of soundings 
taken at the times shown, Figure 3c is a graph of the mixed 
layer depth h and the surface heat flux Q as functions of 
time. The flux is represented by a simple cosine function 
with a zero phase angle. The Monin-Obukov length h is 
calculated as in Eq, (15). 

The first test performed on the TYPE I model was a 
comparison and significance test on the effects of diffusion. 
The regular diffusion coefficient k=l,5 cm 2 sec 1 was consistent 
with deep ocean dynamics as explained by Munk (1966). The 
experiment was to see if the variation of k would change the 
annual temperature distribution, k was arbitrarily raised 
by an order of magnitude over Munk's value and then lowered 
by an order of magnitude. For large diffusion, a rapid 
conduction of heat downward combined with the boundary 
condition of zero heat flux at the bottom would produce a 
tendency toward isothermality , With very small diffusion, 
stronger thermal stratification was expected. In essence 
small or weak diffusion caused temperature changes to be 
affected only by the convective and wind mixing processes. 

This is readily seen in Fig, 4a by the near horizontal 
sections of the 20°C and 21°C isotherms. The obvious 
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differences of Figs, 4a and 4b pointed out the importance of 
choosing a diffusion coefficient that did not overpower the 
other processes while simultaneously providing representative 
diffusion at depths to. which the other processes did not 
extend. In order to accomplish these properties the value 
of 1.5 cm 2 sec^ 1 was retained. 

The next TYPE I experiment was designed to test the 
effects of the mixing parameterization, Computer runs were 
started with the same initial conditions except that the 
mixing parameterization was omitted. Figures 5a and 5b 
represent the basic TYPE I model without wind mixing, The 
striking similarity between Figs. 3a and 3b and Figs. 5a 
and 5b is readily apparent, Convection and diffusion alone 
performed the bulk of the heat transfer in the basic model. 
Wind mixing did not change the temperature profile exten- 
sively. It was therefore concluded that the basic asymmetry 
between spring and fall seen in the previous figures is due 
to the convective adjustment process and not due to mixing, 
Close inspection revealed that the important differences due 
to mixing occurred in autumn which was expected. In autumn 
the ocean was thermally stratified and most vulnerable to 
large downward eddy fluxes due to mixing, Of course the 
thermal stratification was greatest during summer but the wind 
was weak, causing a very shallow mixed-layer depth. In the 
case with no mixing, only diffusion and convection served to 
redistribute heat below the surface, Hence in summer when 
no convection was occurring in the top most layers, the 
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isotherms all tended to deepen at the same rate until autumn 
when convection finally had cooled the ocean down to that 
level, From this point on, convection dominated the vertical 
transport of heat. In the autumn, the deeper isotherms were 
pushed downward further with wind mixing occurring (Fig, 3a), 
than without mixing (Fig. 5a). This was because in the fall 
when the mixed layer was deepening, the requirement that the 
heat flux at -h be downward necessitated that the layer 
beneath -h be heated. In the winter, when the net surface 
flux was upward, the assumption that the flux at -h be 
downward (or zero) required that all levels above -h to cool 
at nearly the same rate, maintaining the stratification. 

The net result was a delay of the abrupt isothermal condition 
seen in Fig, 5a plus a continued deepening of low-level 
isotherms due to strong mixing, 

Use of the Monin-Obukov length as a justifiable mixed- 
layer depth presupposed a stable fluid, During unstable 
conditions, this length was too small, however it was thought 
that convective adjustment would dominate the redistribution 
of heat in the layer. A comparison of a run in which mixing 
was allowed only in the stable case, (Fig. 6a) and mixing 
regardless of stability (Fig. 3a) was made. The definition 
of stability was made directly from the sign of the net 
surface flux. One can observe (Fig, 6a) a seeming disconti- 
nuity in autumn when mixing was abruptly cut off. Although 
the justification for using the Monin-Obukov length did not 
extend to unstable cases, the ocean temperature was still 
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stably stratified, In addition, with the stronger winds 
associated with winter, one expected mixing to at least the 
depth defined by the Monin-Obukov length. When stratified 
conditions are eroded, convective activity was expected to 
dominate and performed thusly, Therefore the decision was 
made to perform wind mixing throughout the year regardless 
of stability. Figure 3a was thus accepted as the TYPE I 
standard against which comparisons were made. 

An additional comparison was made with the model employ- 
ing a mixed-layer depth derived from mixed-layer theory. 

From Denman (1973), when mixing does not extend to the depth 
of the mixed layer, the mixed layer depth is determined by 
the following diagnostic equation; 



Since the numerator was always positive as a consequence of 
the assumption of zero dissipation, one could see the sign of 
h was determined by the denominator, For this purpose, if 
one neglected S ^ , the remaining terms were proportional to 
surface flux. Since h was required to be positive, the 
corresponding surface flux was required to be downward or 
positive. Thus mixing was applicable only in the stable case, 
Figure 7 represents the mixing employing h derived from 
Eq . (24). It was found that this depth was always shallower 

than that from Monin-Obukov. The results compare with Fig. 6a 
When a numerical comparison was made, the results showed that 



h 




(24) 
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the h from Eq , (24) varied from 30% to 80% of h mQ , However, 

for an h > 50 meters, the depth from Eq, (24) was at least 

65% of h , Since mixing was performed only when surface 
mo 

heat flux was downward, (or in April through September), one 
would expect h to be small due to lack of surface winds, 
which were at a minimum during that period, 

From Eq. (15) one can see that there were virtually no 
restrictions on the magnitude of h, In fact, if Q approached 
zero, h would approach infinity. Therefore a restriction on 
the maximum value of h had to be devised. For finite 
differencing purposes, h had to be above the top of the lowest 
layer, A better physical maximum was that h < w + /f, Various 
runs were made using each of the two restrictions. Since no 
problems were evident in either case, the better physical 
restriction that h < w^/f was retained, 

B. TYPE II BOUNDARY CONDITIONS 

As stated before, the TYPE II boundary condition specified 
the air temperature as a function of the time of the year, 

The predicted sea surface temperature from the model was then 
used with the specified air temperature in determining surface 
heat fluxes. Sensible and latent heat were specified 

«2 < T o - T a (t) > 

where Q ^ was 70 ly day" 1 °k , ~ 1 (Haney, 1971), The air 
temperature, T (t), was specified as a cosine function plus 

fir 

a constant such that the net surface heat flux over the 
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annual cycle was close to zero, Longwave back radiation Q D 
was set equal to a constant and was specified as a cosine 
function. Initial conditions of the model, similar to those 
of TYPE I, were isothermal, Also similar to TYPE I runs, 
each run was integrated through 20 years of time when an 
equilibrium was attained, It was thought that if the air 
temperature was significantly warmer or cooler than the ocean 
initially, then gradually the ocean would adjust and even- 
tually reach an equilibrium state. The initial isothermal 
conditions was designed to minimize the transient stage, 

TYPE II experiments were of two general categories; 
those dealing with the extrapolation which defines the ocean 
surface temperature and those dealing with comparisons against 
observations and TYPE I runs, The TYPE II conditions were 
purposely run in order to approximate more closely the 
conditions used in the three-dimensional baroclinic model 
(see section IIB) so the results would reflect the results 
of the large model when mixing scheme was introduced, 

The model predicted the temperature at all ten levels, 
with the 10 meter temperature being the shallowest, This 
predicted temperature T^ represents the average temperature 
when integrated over the depth of the top layer, 20 meters. 
Thus the actual surface temperature could be somewhat warmer 
than T^ , The temperature at 30 m was the predicted value T^. 

If an isothermal layer was assumed to exist to an arbitrary 
depth h' >0 (see Fig, 8) one could derive the following 
expression for the surface temperature; 
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( 25 ) 



rp _ rp (60 + h' ) _ „, ,20 - h\ ' 1 

sfc 1 2 2 ^ 2 ; (h’ + 20) 



where the temperature at 20 meters was 



T 20m 2 



It can be seen that if h'=0, such that no isothermal layer 



existed , 


then 
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Figures 


9, 10, and 


11 


were 


the results from the 


experiments 



corresponding to the extrapolations defined in Eqs, (26a, 26b, 
and 26c) respectively. No extrapolation exhibited undue 
sensitivity even though Fig, 11 displayed a secondary sea- 
surface temperature maximum in the fall, In fact each of 
the above extrapolations displayed this secondary maximum 
but it was not evident in all of the figures due to the 
contouring routine which drew only integer isotherms. The 
magnitude from relative minima to relative maximum in the 
secondary oscillation was 0,3°K, It appears that the 
extrapolations did cause the equilibrium heat content over 
the annual cycle to differ, In Eq . (25), the surface 
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temperature decreases with increasing h , Thus for the 



coldest surface temperature, T ^ =T^ , the ocean exhibited its 
largest heat content over the annual cycle, Remembering that 
air temperature was specified, the warmest surface extrap- 
olation would equilibrate with cooler temperatures at levels 
one and two than the other extrapolations, Thus, the total 
heat content over the annual cycle seemed to be the only 
real difference of the three extrapolations. The value of 
T s £ c in Eq. (26b) was used in the three-dimensional 
baroclinic model, 

Figure 12 shows the case with no mixing using TYPE II 
conditions. The result should be compared with Fig, 11 
since h'=20 meters. The effects of mixing in the TYPE II 
model were significant, Not only did the mixing produce the 
same effects here as in TYPE I models, but it produced the 
secondary maximum. This occurred because in the early fall, 
mixing transported the heat down faster than it was coming 
in at the surface. This caused a short period in early 
September in which the flux at the bottom of the level exceeded 
the surface flux, causing a temporary cooling near the end of 
the warming season. This produced the local minimum followed 



by the secondary maximum in the early fall. The secondary 
maximum occurred because the downward flux into the upper 
level was greater than the downward flux to the second level. 
The downward flux out of the top level apparently decreased 
due to the increasing h. Comparisons between Figs, 11 and 
12 readily show the gradual deepening of the isotherms 
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during summer due to diffusion, However, the local sea- 
surface temperature minimum in fall is not present in Fig, 12 
nor is the delay of isothermal conditions. The conclusion 
was that both these effects were due to the mixing parameter- 
ization , 

The TYPE I conditions (compare Fig, 3a and Fig, 11) 
displayed all the effects due to mixing stated above except 
for the local minimum in fall, The TYPE II boundary condition 
portrayed a much more complex interaction between ocean and 
atmosphere, whereas the TYPE I ocean was one that simply 
reacted to the atmosphere. One might expect that a stratified 
ocean subject to vigorous mixing due to autumn storms would 
cool the upper layers of the ocean through the downward trans- 
port of heat inherent in wind mixing. If the downward surface 
heat flux continued, warming the top layer, it would be 
reasonable to expect the local minimum, This premise 
presupposes strong mixing at a time when air temperature was 
not cooling at a fast rate. The phase relationship between 
mixing (seen in the magnitude of h) and surface flux would 
be critical to the presence of the secondary maximum, 

The secondary maximum was not a phenomena found solely in 
the models shown here, Dorman ejt, al, (1974) found obser- 
vations at Ocean Station "N" to have a late summer local 
minimum and a secondary maximum (Fig, 12). Although the time 
of occurence was not exactly the same in the model and in the 
observation, it must be remembered that the representation 
of air temperature (cosine function) specified in the model 



35 



had no phase angle nor did the wind .from which h was 
calculated. Neither the initial conditions nor the boundary 
conditions were designed to represent Ocean Station "N," but 
the trends in the isotherms were expected to be similar. As 
stated before, the lack of steady state vertical motions 
caused the failure of a permanent thermocline to develop, 

In Fig. 13 , isotherms of 17°C and below could be considered 
part of the permanent thermocline, 

Although in some ways Fig, 10 and Fig, 3 should 
resemble the observations in Fig, 13 in no way were they 
designed to represent that particular ocean station under 
the same conditions. The one^dimensional models shown here 
were for testing only. They were simplified grossly so that 
the complexities of the results from the three dimensional 
model could be more readily understood. 
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IV. THREE-DIMENSIONAL BAROCLINIC MODEL RESULTS 



The three-dimensional baroclinic model began from rest 
with temperature a function of depth and latitude only. The 
model was integrated forward in time 5 years with constant 
forcing; that is using only the annual mean values of 
atmospheric parameters discussed in section lib. This was 
done in an effort to reduce the initialization shocks and 
to "spin up" the currents to a quasi-steady state. Time 
dependency was introduced to the atmospheric parameters on 
the sixth year of integration and was continued for the rest 
of the integration. The mixing parameterization was introduced 
on the ninth year of integration, and the results shown below 
are for the tenth year. 

The data of two points of geographical interest were 
saved every 5 days for later analyses, These were near Ocean 
Station November (30°N, 140°W) and Ocean Station Papa (50°N, 
145°W). Both of these stations have been the subject of 
bathythermograph data collection since the late 1940' s, As 
stated before, Fig, 13 represents the average data collected 
at Ocean Station November over a seven-year interval, The 
baroclinic model in its tenth year of integration produced 
Fig, 14 which shows the annual fluctuation of temperature 
from soundings taken every 5 days. The figure shows several 
aspects encountered in one dimensional test runs, The 
basic asymmetry of the isotherms is present as well as the 
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characteristic deepening of isotherms in summer months due 
to diffusion, However in fall the temperature pattern is 
typical of that seen in the no mixing case, Fig, 12a, The 
absence of a secondary maximum of any magnitude reinforces 
the belief that the mixing process in this model Was very 
weak or completely absent, The conclusion drawn here was 
that the climatological winds provided too shallow a Monin- 
Obukov depth to produce any noticeable effect. 

The permanent thermocline in the model, when compared to 
Fig. 13 , is too diffuse and as a result, the temperatures 
at depth (150m to 280m) were therefore too high. One possi- 
bility is that the vertical diffusion coefficient was too 
large or the vertical velocity "w" was of insufficient magni- 
tude to support a permanent thermocline with a realistic 
gradient. Another possibility is that this deficiency, 
which also exists in other general circulation models 
(Bryan, 1975), may be due to the absence of a vertical eddy 
heat transport by baroclinic meso-scale eddies which are known 
to populate the real oceans, Finally, the magnitude of the 
surface temperatures was about 3°C too large. This may have 
been simply due to the temperature forcing which was the 
surface air temperature averaged over the entire latitudinal 
belt. Clearly this type of averaging would tend to warm 
those normally cool points and cool those normally warm points. 
More precise forcing (individual data at each grid point) 
was unfeasible due to insufficiently accurate data, and 
computer storage restrictions, 
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The maximum and minimum surface temperatures in the 
model's "N" point occurred abdut the first of September and 
the first of March respectively, The observations from 
Dorman could not be compared for the maximum point due to 
the presence of the double maxima. However the middle of 
September seemed plausible if the local minimum had been 
absent. The minimum temperature occurred in the first half 
of April, Thus, although the time duration of the heating 
cycle was approximately correct, the phase lag with respect 
to observations was somewhat small, These discrepancies 
were not considered serious. 

Another point was sampled for its temperature cycle. 

Ocean Station Papa manned by Canadian vessels was selected 
for comparisons. Figure 15 was produced with data from 
the tenth year of integration at the model's point corre- 
sponding to Papa's location. 

Although a time dependent plot similar to Dorman's was 
not available, monthly data was available from Scripps 
Institute of Oceanography in the NORPAX project. The absence 
of a permanent thermocline at this latitude is in accordance 
with observations but the model's temperatures were also too 
warm here. Figure 15 also shows a lack of effects due to 
mixing, Since a time dependent plot of observations was not 
available, phase lag and duration were impossible to compare 
accurately with available data, 
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The three dimensional model, unlike the one dimensional 
models discussed above, had horizontal advection as an active 
heat transport mechanism, The seasonal variation of temper- 
ature advection is likely to be important for NORPAX goals, 

It did not appear to have an important role (compare Fig, 

14 and Fig, 10) but has not been analyzed in this paper. 

This analysis is left for further study, 
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V, CONCLUSIONS 



Although the baroclinic model did not perform up to 
expectations, most causes were clear and correctable. With 
proper initialization, a parameterization of storms which 
would give higher winds, plus an analysis of the diffuse 
permanent thermocline, it was thought that the baroclinic 
model could be used as a tool to study the evolution of the 
temperature anomalies in the North Pacific. 

It has been shown that the climatological winds produced 
a forcing that made deep convective overturning dominant 
over other heat transport mechanisms, Thus on a model using 
climatological atmospheric parameters, the mixing parameter- 
ization may well be ignored or omitted, The model discussed 
in this work was one that could be coupled with an atmospheric 
model in which actual data could be used, thus providing the 
mixing parameterization the needed winds to perform adequately. 

It was thought that the mixing parameterization, through 
the results shown in Section Illb, produced effects on the 
temperature profile characteristic of actual mixing. Thus 
the model fulfilled in a large degree the initial purpose 
of better describing the upper ocean dynamics by incorporating 
the effects of surface-generated wind and convective mixing. 
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Figure 2. Representation of the surface flux boundary condition in the 
TYPE I model. Fluxes are in Langleys per day, 




Figure 3a. The annual cycle of temperature distribution 
using TYPE I boundary conditions and mixing. 
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Figure 3b. Selected soundings from Fig. (3a). 




Figure 3c, Mixed layer depth h in hundreds_of meters and 
surface heat flux Q in cm°k sec ^ x 10~ , 
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Figure 4a, The annual cycle of temperature distribution 

using TYPE I boundary condition with k= , 15cm^sec - ^ . 




Figure 4b. The annual cycle of temperature distribution ^ ^ 

using TYPE I boundary conditions with k =15, Ocm sec - . 
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Figure 5a, The annual cycle of temperature distribution 

using TYPE I boundary condition with no mixing, 




Figure 5b. Selected soundings from Fig, 5a. 
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Figure 6a. The annual cycle of temperature distribution using 
TYPE I boundary conditions with mixing only in 
the Stable Case. 
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Figure 6a. Selected soundings from Fig. 6a. 
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Figure 7. The annual cycle of temperature distribution using 
TYPE I boundary conditions with h calculated from 
Eq. (24), 
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Figure 8. Representation of the extrapolation of sea 

surface temperature using an isothermal depth h', 
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Figure 9. 



Figure 10, 




The annual cycle of temperature distribution using 
TYPE II boundary conditions with T gfc =3/2 T-^-jT^ 




The annual cycle of temperature distribution using 

TYPE II boundary conditions with T _ =7/6T 1 -1/6T Q . 
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Figure 11, The annual cycle of temperature distribution using 

TYPE II boundary conditions with T , =T- . 

s l c 1 
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Figure 12. The annual cycle of temperature distribution 
using TYPE II boundary conditions with no 
mixing . 
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Figure 13, The annual cycle of temperature distribution 

at Ocean Station "N" from Dorman et al (1974), 




Figure 14, The annual cycle of temperature distribution 

predicted by the model from Ocean Station "N. M 
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Figure 15, The annual cycle of temperature distribution 
predicted by the model for Ocean Station "P , " 
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